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We  consider  nonlinear  systems  dx/dt  =  f(x(t ))  where  Df(x(t ))  is  known  to  lie  in  the 
convex  hull  of  L  matrices  Ai,  .  .  . ,  Al  G  RnXU  For  such  systems,  quadratic  Lyapunov 
functions  can  be  determined  using  convex  programming  techniques  [1],  This  paper  de¬ 
scribes  an  algorithm  that  either  finds  a  quadratic  Lyapunov  function  or  terminates  with 
a  proof  that  no  quadratic  Lyapunov  function  exists.  The  algorithm  is  an  interior-point 
method  based  on  the  theory  developed  by  Nesterov  and  Nemirovsky  [2]. 

1.  AN  EQUIVALENT  OPTIMIZATION  PROBLEM 

Quadratic  Lyapunov  functions  of  the  form  xT Px  can  be  determined  by  solving  a  set  of 
matrix  inequalities  for  P: 

A\P  +  P Ak  <  0,  k  =  1, . . . ,  L  .  . 

kI<P<I , 

where  /  is  the  n  X  n  identity  matrix,  and  the  notation  X  >  0  means  that  X  is  positive 
semidebnite.  The  lower  and  upper  bounds  on  P  are  added  for  numerical  reasons.  These 
bounds  limit  the  condition  number  of  P  to  1/k,  which  is  not  very  restrictive  as  long  as  k 
is  a  small  number.  The  upper  bound  also  guarantees  that  the  set  of  feasible  matrices  P 
is  bounded. 

Problem  (1)  can  be  converted  into  an  optimization  problem  by  adding  an  artificial 
variable  t: 

minimize  t 

such  that  —A^P  —  PAk  +  tFk  >  0,  k  =  1, . . . ,  L 

—P  +  tiA+i  +  /  >  0  (2) 

P  +  tFh- 1-2  —  kI  P  0 
t  >  0 

The  optimization  algorithm  requires  a  feasible  solution  to  start  with.  Assume  one  has 
a  symmetric  matrix  P0  as  initial  guess  for  the  matrix  P.  One  can  then  choose  F \  = 
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/  +  A^P0  -f  PoAk ,  k  =  1, . . . ,  L,  Fl+i  =  Po  and  Fl+2  =  (1  +  k)I  —  P0-  This  implies  that 
t  =  1,  P  =  Po  satisfies  the  constraints  (2).  Starting  with  these  initial  values,  one  minimizes 
t.  If  the  minimum  value  of  t  is  zero,  a  solution  to  the  original  feasibility  problem  (1)  has 
been  found.  If  on  the  other  hand  the  minimum  value  of  t  is  greater  than  zero,  we  have 
actually  proven  that  the  original  problem  is  infeasible. 

2.  POTENTIAL  FUNCTION 

A  number  of  very  efficient  methods  have  been  recently  developed  for  solving  optimiza¬ 
tion  problems  over  the  cone  of  positive  definite  matrices.  These  so-called  interior-point 
methods  were  first  introduced  for  linear  programming  by  Karmarkar  [3].  Since  then,  they 
have  been  extended  to  general  convex  problems.  The  first  unified  treatment  can  be  found 
in  [2],  It  is  now  possible  to  solve  problems  involving  large  matrix  inequalities  in  a  rea¬ 
sonable  amount  of  time.  The  worst-case  effort  can  be  shown  to  grow  polynomially  with 
problem  size;  the  performance  is  even  better  in  practice. 

The  most  efficient  interior-point  methods  are  based  on  a  potential  function.  For  prob¬ 
lem  (2)  one  can  use 

L+ 2 

=  glogt  -  logdetAh  (3) 

k= 1 

where 

Xk  =  —A^P  —  PAk  +  tFk,  k  =  l,...L 
Xl+i  =  —P  +  + 

Xl+ 2  =  P  +  tFL+2  —  Kl , 

and  q  is  a  positive  scalar.  The  potential  function  ip(P,t)  consists  of  two  parts.  The  first 
term,  q  log  t,  is  concave  and  rewards  a  decrease  in  t.  The  second  term,  —  J2k=i  log  det  Xk} 
acts  as  a  barrier  for  the  feasible  set. 

Although  the  function  itself  is  not  convex,  it  can  be  shown  that  exp  is  convex  if 
q  >  n(L  +  2).  This  property  has  the  important  consequence  that  local  minima  of  will 
necessarily  be  unique  and  global. 

3.  MINIMIZATION  OF  THE  POTENTIAL  FUNCTION 

Interior-point  methods  use  a  damped  Newton  method  for  minimizing  the  potential 
function  (see  [2]).  Skipping  details,  this  amounts  to  repeatedly  computing  a  step  SP}  St 
by  solving  the  overdetermined  system: 

I  «  X~1/2  ( -A\SP  -  SPAk  +  6tFk)  X“1/2,  k  =  1, . . . ,  L 

/  «  X^+{2(-SP  +  StFL+1)X^+{2 

I  «  X£l{2(SP  +  StFL+2)X£l{2 
ss  t~1St 


-q 


(4) 
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in  a  least-squares  sense.  One  then  takes  a  suitably  scaled  step  in  this  direction. 

If  the  problem  is  feasible,  it  can  be  shown  that  this  Newton  step  reduces  the  potential 
function  at  least  by  an  absolute  constant  in  each  iteration.  This  property  is  crucial  in 
establishing  a  proof  of  polynomiality  of  the  algorithm  [2]. 

4.  DETECTION  OF  INFEASIBILITY 

If  the  potential  function  is  not  unbounded  below,  but  reaches  its  minimum  for  some 
(P*,P)  with  t*  >  0,  we  may  conclude  that  the  problem  (1)  is  infeasible.  A  formal  proof 
of  infeasibility  can  be  derived  from  duality  theory  [2].  It  consists  in  generating  a  set  of 
n  X  n  matrices  Z k  >0,  &  =  1, .  .  .  ,  P  +  2  that  satisfy: 

Tr  Zl+ i  —  kTV  Zl+ 2  <  0 

L 

—  22  +  AkZ2j  —  Zl+ i  +  Zl- |_2  =  0 

k= i 

L+ 2 

^  FkZk  =  !• 

k= 1 

A  set  of  matrices  Z k  satisfying  these  conditions  can  be  derived  from  the  minimizers  (P*,  t*) 
at  a  negligible  computational  cost. 

5.  NUMERICAL  EXAMPLE 

Interior-point  methods  typically  require  a  small  number  of  iterations,  but  each  iteration 
involves  solving  a  large  least-squares  problem.  There  is  no  need,  however,  to  solve  these 
least-squares  problems  exactly.  Approximate  solutions  can  already  yield  descent  directions 
that  are  sufficient  for  polynomiality  [4],  The  use  of  iterative  methods  such  as  conjugate 
gradients  and  LSQR  [5]  to  approximately  solve  the  least-squares  problems  can  therefore 
lead  to  very  considerable  savings  in  computer  time. 

The  numerical  data  for  Fig.  1  were  obtained  as  follows.  We  randomly  generate  20x20 
matrices  SAk,  k  =  1,...,50,  and  construct  problems  of  the  form  (1)  by  taking  Ak  = 
— I  a  6  Ak,  k  =  1,...,P. 

The  problem  will  clearly  be  feasible  for  small  values  of  a  and  infeasible  for  larger 
values.  Figure  1  shows  the  total  cpu-time  (on  a  SUN  4/670)  and  the  number  of  iterations 
for  different  values  of  a.  From  these  data,  we  note  that  the  algorithm  is  fast  if  the  problem 
is  clearly  feasible  or  infeasible,  and  becomes  slower  as  we  approach  the  boundary  between 
feasibility  and  infeasibility. 

In  this  example,  the  optimization  problem  (2)  has  211  unknowns  (the  scalar  t  plus  the 
elements  of  the  symmetric  matrix  P).  The  least-squares  problem  (4)  that  has  to  be  solved 
in  each  iteration  has  size  10921  X  211. 
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cpu-time  (seconds)  no.  of  outer  iterations 


Figurt  1.  Total  cpu-time  and  number  of  iterations  for  the  experiment  described  in  Sec¬ 
tion  5.  The  instances  marked  with  ’ x’  were  feasible;  the  instances  marked  with  ’o’  were 
infeasible. 
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